function [xst, info] = NMBB(fun,x0)
% Barilai-Borwein descent algorithm.
    eps = 1.0e-6 ; 
    iterM = 1000;
    tol = 1.0; 
    iter = 0; 
    ls = 0; 
    f_num = 0; 
    g_num = 0;
    m = 3;  %nonmonotone ls parameter
    start_time=tic();
while (tol>=eps) && (iter<iterM)
     if (iter == 0)
        g0 = grad(fun,x0);
        g_num= g_num+1;
        xk =  linesearch(fun,x0,-g0, g0); 
        ls = ls+1;
    end
    gk = grad(fun,xk);
    g_num = g_num+ 1;
    tol = sqrt(gk'*gk);     
    sk = xk - x0; 
    yk = gk-g0;
    alpha_bb = (sk'*sk)/(sk'*yk);
    x0 = xk;  
    xt = x0 - alpha_bb*gk;
    g0 = gk;
    ft = feval(fun,xt); 
    f_num = f_num +1;
    fval = zeros(m-1,1);     
    for i = 1:(m-1)
       if (mod(iter,m)== i)
           fval(i) = ft;
        end
    end
    fmax = max(fval);
    if (mod(iter,m)==0)
        [xt,ft] = nmlinesearch(fun,x0,-gk,gk,fmax);
        ls = ls +1;
    end
    xk = xt; 
    fk = ft;
    iter = iter +1; 
end
        cput = toc(start_time);  
        
        xst = xk;
        info.iter = iter;
        info.time = cput;  
        info.f_num = f_num;
        info.g_num = g_num;
        info.f = fk;
        info.g = tol;  
        info.ls = ls;
        info.bb = iter-ls;
        info.rator = (iter-ls)/iter;
end